Ground-State Phase Diagram of the Two-Dimensional Extended Bose-Hubbard Model 
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We investigate the ground-state phase diagram of the soft-core Bose-Hubbard model with the 
nearest-neighbor repulsion on a square lattice by using an unbiased quantum Monte Carlo method. 
In contrast to the previous study[P. Sengupta et. al, Phys. Rev. Lett. 94, 207202 (2005)], we 
present the ground-state phase diagrams up to large hopping parameters. As a result, in addition 
to the known supersolid above half-filling, we find supersolid even below and at half-filling for 
large hopping parameters. Furthermore, for the strong nearest-neighbor repulsion, we show that 
the supersolid phase occupies a remarkably broad region in the phase diagram. The results are 
in qualitative agreement with that obtained by the Gutzwiller mean-field approximation [M. Iskin, 
Phys. Rev. A 83, 051606(R) (2011) and T. Kimura, Phys. Rev. A 84, 063630 (2011)]. 
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I. INTRODUCTION 

Supersolid has attracted great interest for a long time 
as a fascinating quantum state that has superfluidity and 
solidity simultaneously. In the early theoretical works 
by Andreev and LifshitzpJ, and by Chester 0, they pro- 
posed a scenario that supersolid might appear when zero- 
point defects in solid such as 4 He undergo Bosc-Einstein 
condensate at low temperatures without destroying the 
crystal structure. After several decades, a discovery was 
made in 2004 by Kim and Chang, [|. In their exper- 
iments on solid 4 He, they observed nonclassical rota- 
tional inertia associated with superfluidity in the solid. 
After the discovery, further theoretical or experimental 
works 0-0] provided the evidence that it is different from 
a bulk supersolid of the Andreev-Lifshits-Chester sce- 
nario. The superfluidity in solid 4 He seems to appear due 
to the extended defects such as grain-boundaries d, Q or 
dislocations [loj. 

In contrast to the supersolid in continuous spaces, su- 
persolid in lattice systems has been a promising candi- 
date recently. This is based on the recent experimental 
development of optical lattice systems (liTtl4l ] . Ultra-cold 
Bose gases trapped in optical lattice are ideal systems to 
realize the Bose-Hubbard models [l5|. From the inten- 
sive theoretical and numerical studies, the existence of 
supersolid phases has been established in the extended 
Bose-Hubbard models ^6H3l1 |. Most of the supersolids in 
lattice systems are achieved by doping particles or holes 
into insulating solid states at commensurate filling fac- 
tors. If doped defects delocalize and Bose-Einstein con- 
densate against a phase separation, supersolid appears by 
the the Andrccv-Lifshits-Chcstcr scenario. Thus, result- 
ing supersolids are stabilized at incommensurate filling 
factors. 

One of the simplest models to study supersolids is the 
soft-core Bose-Hubbard model with nearest-neighbor re- 
pulsions. By accurate quantum Monte Carlo calculations 
on this model, chcchcrboard supersolid phases have been 



found on a ID chain[22j, a 2D square lattice[21(, and 
a 3D simple cubic latticeJH [H In the ID and 

2D cases, supersolid regions are found only above half- 
filling (interstitial supersolid). In contrast, in the 3D 
case, supersolids are also found even below and at half- 
filling for large hopping parameters (vacancy supersolid 
and commensurate supersolid respectively) [261 [30L l32j . 
Especially, the presence of supersolid at the commen- 
surate filling factor 1/2 is fascinating as an exceptional 
supersolid without any doping, although such supersolid 
regions have not been found so far in ID and 2D. There- 
fore, it is a question why there is a discrepancy between 
2D and 3D systems. 

Recent works based on the Gutzwiller mean-field ap- 
proximation have provided some interesting results on 
the ground-state phase diagram of the mo del [H, [34|, in- 
cluding a possible answer to the above question. In the 
ground-state phase diagram presented in Ref. (33|, the 
author found a supersolid phase below and at half-filling. 
Since it was found more clearly for larger hopping param- 
eters, he suggested that the absence of such supersolid re- 
gions in the 2D quantum Monte Carlo study[2l[ might be 
due to the not so large hopping parameter. Since the re- 
gions for the supersolid below and at half-filling are much 
smaller than that above half-filling and the mean-field 
approximations tend to overestimate the region of the 
supersolid phase [2(| HH, the existence of such supersolid 
regions is a subtle problem. As discussed in Ref. [33| . 
more precise treatments are desirable to conclude the ex- 
istence of 2D supersolid phase below and at half-filling, 
because the Gutzwiller approximation becomes more ac- 
curate in higher dimensions and particle densities. 

The other interesting result presented in Ref. [34| is 
on supersolid phases for strong nearest-neighbor repul- 
sion. The ground-state phase diagrams show that, as 
the nearest-neighbor repulsion increases, the supersolid 
phase expands up to large hopping parameters in the 
phase diagram. Especially, the 2D case of this result 
might be most important, because it has the possibil- 
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ity of realizations in quasi-2D dipolar Bose gases whose 
the dipoles are polarized along the z-axis[36j. Therefore, 
from the viewpoint of experiments, we also need to deter- 
mine the more precise phase boundaries in the 2D system 
and check the accuracy of the phase diagram. 

In this paper, motivated by the results of the 
Gutzwiller treatment, we investigate the ground-state 
phase diagram of the extended Bosc-Hubbard model on a 
square lattice by numerically exact quantum Monte Carlo 
simulations. The paper is organized as follows. In Sec. 
ITU we describe the model discussed in this papers and the 
quantum Monte Carlo method we used. Sec. Mil presents 
the ground-state phase diagrams in the grand canonical 
ensembles. These phase diagrams include up to the third 
insulating lobes. Within this region, we confirm that our 
ground-state phase diagrams are qualitative agreement 
with those obtained by the Gutzwiller approximation. In 
Sec. IIV1 we study quantum phase transitions and explain 
the procedure of determining the phase boundaries pre- 
sented in the previous section. In Sec. |V] we investigate 
the supersolid phase at half-filling by obtaining results 
for the canonical ensembles. By showing a ground-state 
phase diagram at half-filling, we confirm that the super- 
solid phase is easily found for large hopping parameters. 
Finally, in Sec. IVI1 we summarize our results. 



II. MODEL AND METHOD 

The model considered in this paper is the soft-core 
Bosc-Hubbard model with nearest-neighbor repulsions on 
a square lattice. The Hamiltonian is given by 

H = -t ]T {b\bj + h.c.) -V^2ni + —^2 n ii n i - l ) 

(id) 1 1 

+Vj2 n * n J- (!) 

(id) 

Here, b\{b/) is the bosonic creation (annihilation) oper- 
ator on site i, and n, is the particle number operator 
defined as rii = b\bi. The summation is taken over 
all pairs of nearest-neighbor sites. For a square lattice, 
the coordination number z equals 4. Furthermore, t is 
the hopping parameter, p is the the chemical potential, 
U is the on-site interaction, and V is the nearest-neighbor 
interaction. In this paper, we consider the case where the 
interactions are repulsive (U, V > 0). In our simulations, 
we treat N = L x L systems with the periodic boundary 
condition. 

In the classical limit t/U = 0, the ground-states are 
known and simple [2ll l30l . . When the nearest-neighbor 
repulsion satisfies zV/U < 1, the ground states are 
checkerboard solids at filling factors p = 1/2, 3/2,..., and 
uniform Mott-insulators at p = 1, 2, ... . To characterize 
each state, we can label it as (n^,ns) which represents 
a pair of particle numbers on the two sublattices A and 
B. Without loss of generality, we assume that Ha > n>B- 



Based on this notation, the ground states are labeled as 
(1,0), (1,1), (2,1), (2,2), ... at p = 1/2, 1, 3/2, 2, ... re- 
spectively. In contrast, for zV/U > 1, all ground states 
are checkerboard solids. The states are labeled as (1,0), 
(2,0), (3,0), (4,0), ... at p = 1/2, 1, 3/2, 2, ... respec- 
tively, and the transition from p = n/2 to (n + l)/2 takes 
place at {p,/U) c = n, when the chemical potential is in- 
creased. Therefore, zV/U — 1 is a critical point for p > 1 
in the classical limit. When the finite t/U is introduced, 
the critical point (zV/U) c — 1 is shifted to slightly larger 
values due to quantum fluctuation. 

To investigate the properties of the model for finite val- 
ues of t/U , we used an unbiased quantum Monte Carlo 
method. The formulation we used is based on the Feyn- 
man path integral representation. In the representa- 
tion, the (i-dimensional quantum system is mapped to 
the (d+ l)-dimensional classical systems. In the mapped 
systems, each configuration is called world-line with d- 
dimensional space axises and one-dimensional imaginary 
time axis. Based on this representation, we sample the 
world-line configurations according to the Markov chain 
Monte Carlo. To update the configurations, we used the 
worm- type algorithm [3tM4Q| ] . 



III. GROUND-STATE PHASE DIAGRAM IN 
THE GRAND-CANONICAL ENSEMBLE 

In this section, we present ground-state phase diagrams 
in the zt/U-p/U plane. The recent Gutzwiller mean-field 
study suggested that the supersolid phase might exist 
even below half-filling for large hopping parameters [33j . 
In addition, the other work provided the results that the 
ground-state phase diagram have qualitatively different 
structures between weak nearest- neighbor repulsions and 
strong nearest- neighbor repulsions [34j. Remarkably, in 
the latter case, the supersolid phase seems to occupy 
very large region in the phase diagram. To confirm these 
results by numerically exact quantum Monte Carlo cal- 
culations, we show the ground-state phase diagrams at 
zV/U = 1 and zV/U = 1.5 in Sees. HIS and |TiTB] re- 
spectively. 



A. Ground-state phase diagram at zV/U = 1 

In Fig. Q] (a) , we show the ground-state phase diagram 
at zV/U = 1 in the zt/U-p/U plane. To detect each 
phase, we measured the particle density p — 1 / 'N $^(7it) i 
the superfluid stiffness p s = (W 2 )/(2dt/3L d - 2 ), and 
the structure factor S(k) = l/N 2 Y,i,je ik - ri '((n i nj) - 
(rii) 2 ). Here, (■ • • ) is the thermal average, and W de- 
notes the winding number vector in the path integral 
representation [4lj. (3 represents the inverse temperature 
defined by (3 = 1/T, d is the dimensionality of system 
that is equal to 2 in this paper, k is the wave vector, and 
Tij indicates the relative position vector between sites i 
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FIG. 1: (Color online) (a) Ground-state phase diagram of the 
extended Bose-Hubbard model on a square lattice at zV/U = 
1. Red circles indicate boundaries for the insulating lobes. 
Blue squares represent the SS-SF boundary. The inset is the 
enlarged view of the region around the tip of the first CB lobe. 
Error bars are drawn but most of them are much smaller than 
the symbol size (here and in the following figures). Black line 
is the boundary between the empty region and the SF that 
can be obtained analytically. Other lines are used to guide 
the eyes, (b), (c), and (b) Schematic configurations for the 
insulators at p — 3/2, p = 1, and p — 1/2 respectively. Each 
red circle represents one particle on the sites. 




0.49 
0.48 
0.47 



zt/U = 0.328 




zV/U=l 






T/t. = 0.05 








' SS 


■ □ L = 48 






♦ O £-64 






* A L = 80 






• o L = 96 


I'll! 


1 


"i| 


ft E 


H 

S(ir, it) : 



0.32 0.33 



li/U 



FIG. 2: (Color online) (a) and (b) Physical quantities as 
functions of p/U at (zt/U, zV/U, T/t) = (0.12,1,0.05) and 
(zt/U, zV/U, T/t) = (0.328,1,0.05) respectively. Shaded re- 
gions indicate the supersolid state where p s and S(n, it) take 
finite values simultaneously. 



and j. In our phase diagram up to p/U < 3, in ad- 
dition to a conventional supcrfmid phase(SF), there are 
three insulating lobes at p = 1/2, p = 1, and p = 3/2. 
Schematic configurations are shown in the Fig. [1] (b), 
(c), and (d), respectively. The lobe at p = 1 is a uniform 
Mott-insulating phase (MI), and the others at p = 1/2 



and p = 3/2 are checkerboard- type solid phases (CB) 
characterized by finite value of S(ir, tt). We also con- 
firm the presence of supersolid phases (SS) around the 
insulating CB lobes. The determinations of the phase 
boundaries are explained in detail in Sec. IIVI 

To show the existence of each phase, we plot 
p/U dependence of the measured quantities at 
(zt/U,zV/U,T/t) = (0.12,1,0.05) and (0.328, 1, 0.05) in 
Fig. [2] (a) and (b) respectively. In the case of the small 
hopping parameter zt/U = 0.12 in Fig. H](a), SS phases 
exit above p = 1/2 and around p = 3/2. When parti- 
cles are removed from the checkerboard solid at p = 1/2, 
possible supersolid is unstable against a phase separa- 
tion as known by strong-coupling argument (2lJ . In con- 
trast, for the larger hopping parameter zt/U = 0.328 in 
Fig. [5] (b), we find that SS phase are present even be- 
low half- filling. As seen in the inset of Fig. [1] (a), the 
SS phase covers the tip of the first CB lobe. This re- 
sult suggests that the supersolid can be also stabilized 
at half-filling. In Sec. [V] we present direct evidence 
for supersolid at half-filling by obtaining results for the 
canonical ensemble and excluding possible phase separa- 
tions. In addition to the SS around p = 1/2, the other 
SS phase around p = 3/2 more clearly covers the tip of 
the corresponding insulating CB lobe. Therefore, the su- 
persolid seems to be stabilized even at p = 3/2. The 
present 2D ground-state phase diagram is in qualitative 



and that obtained by the 
HI, HH]. However, we find 



agreement with that in 3 DI32I 
Gutzwiller approximation |l 6l 
that the supersolid regions clearly become smaller as the 
dimensionality decreases. 



B. Ground-state phase diagram at zV/U = 1.5 

For strong nearest-neighbor repulsions, all insulat- 
ing states are checkerboard solid states and, thus, the 
ground-state phase diagram are quite different from that 
for weak nearest-neighbor repulsions. In Fig. [3] (a), we 
present the ground-state phase diagram at zV/U = 1.5 in 
the ground-canonical ensemble. In contrast to the phase 
diagram at zV/U = 1, all three insulating Mott lobes 
are actually the checkerboard solid ones. The schematic 
configurations at p = 3/2, 1, and 1/2 are shown in Fig. 
[3] (b) , (c) , and (d) respectively. Compared with the case 
of zV/U = 1, the insulating lobes extend up to larger 
hopping parameters. This result is reasonable, because 
the strong nearest-neighbor repulsion favors the checker- 
board solid state. The remarkable point is that the con- 
nected SS phase exits, surrounding all the CB lobes. The 
SS phase occupies a broad region up to large hopping pa- 
rameters, and the phase boundary behaves linearly. Our 
result is still in qualitative agreement with that obtained 
by the Gutzwiller approximation (34j. However, the su- 
persolid region is apparently smaller. 

To support the results, we plot the measured quantities 
as functions of p/U at (zt/U, zV/U, T/t) = (0.2, 1.5, 0.05) 
and (0.6, 0.15, 0.05) in Fig. @] (a) and (b) respectively. 
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FIG. 3: (Color online) Ground-state phase diagram in the 
zt/U — fi/U plane at zV/U = 1.5. (b), (c) and (d) Schematic 
configurations for the insulators at p = 3/2, p = 1, and p — 
1/2 respectively. 
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FIG. 4: (Color online) (a) and (b) Physical quantities as 
functions of p/U at (zt/U, zV/U, T/t) = (0.2,1.5,0.05) and 
(zt/U,zV/U,T/t) = (0.6,1.5,0.05) respectively. 



need different treatments to determine the phase bound- 
aries. In the following three subsections, we explain the 
treatments for each phase boundary. 



In Fig. [4] (a), there are three plateaus atp=l/2,l,3/2, 
where S(tt,7t) takes finite value. These plateaus corre- 
spond to the CB phases. Between these regions, S(ir, 7r) 
and p s take finite value simultaneously, indicating the 
SS phase. In contrast, just below p = 1/2, there is no 
SS phase and we observed a clear discontinuity in the 
particle density again. Just below p = l,and 3/2, the 
slopes in the particle density are very steep. However, 
compared with that below p = 1/2, possible disconti- 
nuities are not so clear. Thus, the CB-SS transitions 
might be wcakly-first-order or second-order at this pa- 
rameter. When the hopping parameter becomes smaller, 
we confirmed that the slopes become steeper, suggest- 
ing the presence of a first-order transition predicted by 
the strong coupling arguments 2l| . For larger hopping 
parameter as in Fig. @] (b), all the insulating plateaus 
disappear. In contrast, the SS phases are connected and 
occupy all the region for large chemical potentials. 



IV. QUANTUM PHASE TRANSITIONS 

In this section, we study quantum phase transitions 
and explain how the phase boundaries are determined. 
There are three different kinds of quantum phase tran- 
sitions in terms of symmetry breaking: the transition 
between two phases with different broken symmetries 
(the CB-SF transition), the superfluid transition that in- 
volves the gauge symmetry (the CB-SS transition and the 
MI-SF transition), and the checkerboard-order transition 
where the translational symmetry is breaking (the SS-SF 
transition). Since these quantum phase transitions have 
different properties related to the broken symmetries, we 



A. Solid-superfluid transition 

We begin with the CB-SF transition that appears at 
the lower boundary of the first CB lobe. As observed 
in Figs. [2] (a) and 0] (a) as well as the previous quan- 
tum Monte Carlo works pli |30|. there are finite jumps in 
the particle density at the boundary, indicating a first- 
order transition. This result can be understood from an 
argument on the broken symmetries in each phase and 
the standard Landau- Ginzburg- Wilson paradigm. In the 
CB phase, the broken symmetry is the Zi associated to 
the broken translational symmetry. On the other hand, 
in the SF phase, the U(l) gauge symmetry is broken at 
zero temperature. (Note that, at finite temperatures in 
two dimensions, the SF phase shows not the long-range 
order, but the quasi-long range order.) According to 
the Landau- Ginzburg- Wilson paradigm, a transition be- 
tween two phases with different broken symmetries re- 
sults in a first-order transition or intermediate region 
where both symmetries are broken simultaneously. Since 
an intermediate supersolid phase is absent at the bound- 
ary, the direct CB-SF transition should be a first-order. 
Thus, we simply determined the phase boundary from 
the position of the finite jump in the particle density. 



B. Solid-supersolid transition and 
Mott-insulator-superfluid transition 

At the CB-SS boundaries and MI-SF boundaries, the 
quantum phase transitions are the insulator-supcrfluid 
ones. As for the value of the dynamical critical exponents 
z c , two possibilities are expected: generic transition with 
z c = 2 and special transition with z c = 1 [42j . Because of 
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the difference, we have to determine the transition points 
in different manners. 

The generic transitions are driven by adding/removing 
a particie to/from the insulating phases. In this case, 
the phase boundary can be determined from the finite- 
size scaling analysis of p s for quantum critical points with 
~ c = 2[43||. However, it can also be determined more sim- 
ply from the the zero-momentum Green function G(p — 
0, r) H^] . In the worm algorithm, the zero- momentum 
Green function can be obtained by measuring the Mat- 
subara Green function G(rj,r) = (T T &j(r)&o(0)}- Here, 
T T indicates the time-ordering operator on the imagi- 
nary time r, and 6j(r) is defined by bi(r) = e TH bie~ rH . 
From the asymptotic exponential decay G(p = 0, r) — > 
Z + e~ A + T (t ->■ +00) [Z_e A - r (t -> -00)], we can esti- 
mate the energy gap A + (A_) for creating single parti- 
cle (hole) excitation with p = in the insulating phases. 
In the ground-canonical ensemble, the energy gap corre- 
sponds to the distance between the observed point and 
the phase boundary in the p direction. Thus, we deter- 
mined the phase boundary from the energy gap. Fig. [5] 
shows an example of estimating the energy gap A + in 
the first CB lobe. 




Tt 

FIG. 5: (Color online) Extraction of the energy gap A(+) from 
the zero-momentum Green function G(p — 0, r) in the first 
CB lobe. Solid circles denotes the results obtained by our 
simulation, and the line represents the exponential fit. The 
inset shows the extrapolation of the obtained A (red squares) 
to the thermodynamic limit. 

In contrast to the generic transition, the special transi- 
tion is driven by delocalizing quantum fluctuation. This 
transition occurs at the tip of insulating lobes with fixed 
fx/U. The tip corresponds to a multicritical poi nt where 
z c equals 1 due to a particle-hole symmetry [42[. There- 
fore, to determine the critical point close to the tip in 
the inset of Fig. [TJ we performed the finite-size scal- 
ing analysis of p s for quantum phase transitions with 
z c = 1. In this analysis, the scaling form is given by 
p s ld+z c -2 _ f(SL 1 / u ,/3/L Zd ), where v is the critical ex- 
ponent of the correlation length, S denotes the distance 
from critical points as S = zt/U — (zt/U) c , and / is a 
scaling function. In the present case of d = 2 and z c = 1, 
the value of d + z c ~ 2 equals 1. Therefore, p s L should 
cross at the critical point for different system sizes with 




0.328 0.3285 0.329 0.3295 -0.1 -0.05 0.05 0.1 



zt/U 5L 1 '" 

FIG. 6: (Color online) (a) Plots of p a L as functions of t/U 
near the tip of the first CB lobe. The vertical dashed line 
is placed at the quantum critical point (zt/U) c = 0.32888(8) 
that is estimated from the crossing point, (b) Scaling plots of 
p 3 L. 

fixed fj/L and we can simply estimate it from the cross- 
ing point. Fig. [6] (a) shows one example of this estima- 
tion. In this figure, we estimated the critical point as 
(zt/U) c = 0.32888(8) for p/U = 0.331 that is very close 
to the tip. 

To clarify the universality class of the special transi- 
tion, we proceed to perform the finite-size scaling analysis 
of p s L. In the case of z c = 1, the effective dimension be- 
comes d + z c = 3. Since the breaking symmetry in this 
transition is related to the global U(l) symmetry, this 
quantum phase transition is expected to belong to the 
3D XY universality class. Using the critical exponent 
v = 0.67155 of the 3D XY universality class |46| and 
the dynamical critical exponent z c = 1, we plot p s L as 
a function of bL x l v in Fig. [5] (b). In the figure, we suc- 
cessfully observe the data collapse for large system sizes, 
supporting the validity of the present analysis. 



C. Supersolid-superfluid transition 

Finally, we explain the SS-SF boundaries. The SS-SF 
transition is the checkerboard-solid transition related to 
the Z2 symmetry breaking of the translational symme- 
try. For this quantum phase transition, the critical point 
can be determined from the Binder ratio g defined by 
g = l/2[3— (m 4 )/(m 2 ) 2 ] Here, m indicates the order pa- 
rameter defined by m = l/N^2 i nie l ^' Vi with A; = (ir,ir). 
The scaling form for g is given by g = f{8L 1 / v , j3/L Zc ), 
where 6 = zt/U - (zt/U) c or p/U - (p/U) c . Therefore, 
g for different system sizes should cross at the critical 
point. As a working hypothesis, we assume that the dy- 
namical exponent z c equals 1. In Fig. [3(a), we show the 
p/U dependence of g at zt/U — 0.24 and fit = 0.5i. As 
can be seen in the figure, g actually crosses at a point for 
different system sizes. From the crossing point, we esti- 
mated the quantum critical point as (p/U) c = 0.08455(5) 
for (zt/U,zV/U) = (0.24,1). 

To check the consistency of our analysis and clarify 
the universality class, we analyzed scaling behaviors of 
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FIG. 7: (Color online) (a) Estimation of the SS-SF bound- 
ary from a crossing point of the Binder ratio g for different 
system sizes, (b) and (c) Finite-size scaling plots of g and 
S(ir, 7r)L 2/3c/l/ respectively. 



cal points, as seen in Fig. [T] (a), and, we found that the 
finite jump becomes larger. 
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FIG. 8: (Color online) Finite jump in the particle den- 
sity at the SS-SF boundary for a small hopping parameter 
zt/U = 0.04. Dashed vertical lines are used to separate dif- 
ferent phases. In the classical limit zt/U = 0, the particle 
density changes discontinuously from 1/2 to I at (fj,/U) c = 1- 



S(ir, 7r) as well as g. The scaling form for 5(7r, 7r) is 
given by S{n, n)L 2 ^ u = f{8L 1 / v ,(3/L z '=), where f3 c is 
the critical exponent of the order parameter. Since the 
effective dimension is d + z c = 2 + 1 = 3 and the bro- 
ken symmetry is Z 2 symmetry, the quantum phase tran- 
sition is expected to belong to the 3D Ising universality 
class. Thus, using the critical exponents v = 0.63001 and 
2j3 c /v = 1.03627 of the 3D Ising universality class (47j. we 
plot g and S(ir, 7r)L 2,9c / I/ as functions of 5L 1 /" with fixed 
(3/L in Fig. 0(b) and (c) respectively. As can be seen in 
the figure, the data collapses for large system sizes agrees 
with the expected scaling behavior. 

Exceptional determination of the SS-SF boundaries 
was made for small hopping parameters zt/U < 0.08 
at zV/U = 1, because we observed clear finite jumps 
in the particle density. Fig. [8] shows a jump at the SS- 
SF boundary, indicating a first-order transition. Simi- 
lar discontinuities have been also found in the previous 
quantum Monte Carlo study 21[ . In this region, we de- 



termined the boundary from the position of the jump 
at low temperatures. The discontinuities of the SS-SF 
boundaries seem to be connected to ones of the CB-MI 
boundaries in the classical limit zt/U = where the par- 
ticle density changes discontinuously from 1/2 to 1, 1 
to 3/2,... at the critical points (fj,/U) c — 1,2,... respec- 
tively. In fact, when the hopping parameter is smaller, 
the SS-SF transition points approach the classical criti- 



V. COMMENSURATE SUPERSOLID PHASE 

Most supcrsolids are realized by adding/removing par- 
ticles to/from a commensurate insulating solid. When 
doped defects dclocalizc against phase separations and 
give rise to superfulidity on solid, a supersolid state ap- 
pears. In contrast to this superolid, the situation of the 
supersolid at commensurate filling factors is different, be- 
cause any dopants are absent. In this section, by obtain- 
ing simulation results in the canonical ensemble, we in- 
vestigate supersolid exactly at the commensurate filling 
factor p = 1/2. To obtain results in the canonical ensem- 
ble with the grand-canonical method, we performed the 
following procedures. We first estimated the chemical 
potential that corresponds to the desired particle den- 
sity with high accuracy. Then, we performed simulations 
at the obtained chemical potential and used only samples 
whose particle density is exactly equal to the desired one. 
With this method, in Sec. IV Al we obtain direct evidence 
for supersolid at half-filling, excluding the possibility of 
phase separations. In the following Sec. IV Bl we present 
the ground-state phase diagram at half-filling. The ob- 
tained phase diagram shows that the supersolid phase 
can be found more clearly as the nearest-neighbor repul- 
sion zV/U increases, as the suggestion by the work based 
on Gutzwiller approximation [33j. 
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A. Supersolid at half- filling 

In this subsection, we explicitly show the presence of 
supersolid at half-filling. In Fig. [SJ we plot p s and S(tt, it) 
as functions of the temperature at half-filling. At low 
temperatures, both p s and S(ir, it) have finite values, in- 
dicating a supersolid state. To exclude the possibility 
of phase separations, we show a snapshot of the typical 
configuration in Fig. 1101 In our snapshots, we do not 
find any macroscopic phase separations. Instead, we can 
see that the checkerboard solid has microscopic defects 
(intersitials or vacancies), suggesting the superfluidity is 
caused by dclocalizing defects in the same way as the or- 
dinary supersolids. However, the origin of defects seems 
to be different from the ordinary one, because it is real- 
ized without any change from the commensurate filling 
factor. Since the CB-to-SS transition at half-filling corre- 
sponds to the special transition at the tip of the CB lobe 
in the grand-canonical phase diagram, it is driven not 
by adding or subtracting a particle, but by delocalizing 
quantum fluctuation. Therefore, it is reasonable to inter- 
pret the origin of defects as unbound interstitial-vacancy 
pairs due to the delocalizing quantum fluctuation [f|. 




T/U 

FIG. 9: (Color online) Finite-temperature dependence of p 3 
and S(ir, n) exactly at half-filling. 

Melting of the supersolid occurs through two successive 
finite-temperature transitions, namely supcrfluid transi- 
tion and solid transition. Each critical temperature can 
be determined as follows. We first consider the super- 
fluid transition. In Fig. [§J we can observe the strong sys- 
tem size dependence of p s above the superfluid region, 
which is characteristic of the Kostcrlitz-Thouless(KT) 
transition[48|, To determine the critical temperature 
of the KT transition, we make the x 2 n t to the critical 
form for the squared winding number [HoL l5lj . Specif- 
ically, the squared winding number follows the scaling 
from of {n/4)(W 2 } = 1 + [2hi(L / Lq)]- 1 at the critical 
point. Here, Lq is the only free parameter. For each 
temperature, we make the x 2 fit to the critical form and 
measure the x 2 . Finally, we can obtain a critical temper- 
ature as the temperature that minimizes the value of x 2 - 
The result is shown in Fig. Qj](a). From this analysis, we 
estimated the critical temperature of the KT transition 




FIG. 10: (Color online) Snapshot of supersolid at half-filling. 
This shows a typical configuration in a real space at some 
particular imaginary time. The parameters are chosen at 
(L,zt/U,zV/U,T/U) = (32,0.33, 1,0.008). Each site are do- 
nated as a square. Empty, blue, and red squares indicate 
empty sites, singly-occupied sites, and doubly-occupied sites 
respectively. 



as (T/U) c = 0.0170(5). 

Next, we determined the critical temperature of the 
checkerboard-solid transition from the structure factor. 
For finite-temperature phase transitions, the scaling form 
is given by S(ir, ■k)L^-/ v = f{SL l / u ), where 8 = (T/U)- 
(T/U) c . Since the transition is related to the Z2 symme- 
try breaking, we expect that the critical exponents 2j3 c /v 
and v equal 1 /4 and 1 respectively for the 2D Ising uni- 
versality class. When this is the case, S(tt, ■K)L 2 ^ a l v for 
different system sizes should cross at a critical tempera- 
ture. Fig. [TT](b) shows the result. In the inset, to check 
the consistency on the critical exponents, we present the 
result of the scaling plots that shows the excellent data 
collapse. Therefore, we obtained the critical temperature 
of the checkerboard-solid transition as (T/U) c = 0.066(1) 
from the intersecton of S(tt, ir)L 2 P"/ u . 

(a) (b) 



0.3 




FIG. 11: (Color online) Determinations of the two critical 
temperatures in the supersolid state, (a) Values of \ 2 (solid 
squares) for each temperature. At the critical temperature, 
the value of x 2 is expected to be minimized, (b) Intersection 
of the structure factor S(tt, 7r)L 2/3c/ ' 1 ' for different system sizes. 
The position of the intersection corresponds to the critical 
temperature of the checkerboard-solid transition. In the inset, 
we present the data collapse of the scaling plots. 
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B. Ground-state phase diagram at half-filling 

In the previous quantum Monte Carlo study (2lj. super- 
solid phase has not been found at half-filling for zt/U — 
0.2 [2l|. According to the results from the Gutz wilier ap- 
proximation, this might be because the hopping param- 
eter is not sufficiently large for supersolid to be found 
clearly at half-filling 33] . In this subsection, to confirm 
this suggestion, we clarify the parameter dependence of 
the supersolid region at half-filling. 

In Fig. [I2l we present the ground-state phase dia- 
gram at half-filling in the zt/U-zV/U plane. The phase 
boundaries are determined from the position of an in- 
tersection of g or p s L for different system sizes with the 
assumption that z c equals 1. Fig. [T3] shows a result 
at zV/U = 1. In the figure, we obtained the quantum 
critical points for the CB-SS transition and the SS-SF 
transition as (zt/U) c = 0.32888(8) and 0.33332(8) re- 
spectively. Note that the critical point for the CB-SS 
transition at p = 1/2 agrees with that obtained from the 
grand-canonical cnscmblc(Scc. IIVB[) . In our phase di- 
agram, the supersolid region is much smaller than that 
obtained by the Gutzwiller approximation 33] . However, 
qualitative behaviors of the phase boundaries agree with 
the Gutzwiller results. As the nearest-neighbor repulsion 
zV/U increases, the CB phase expands up to larger hop- 
ping parameters zt/U . The SS phase also extends for 
large nearest-neighbor repulsions and hopping parame- 
ters. In contrast, for the small hopping parameters in- 
cluding zt/U = 0.2, the two phase boundaries are very 
close to each other. Thus, we conclude that the reason 
why the SS phase was not found at half-filling in the pre- 
vious quantum Monte Carlo result [2lJ is that the hopping 
parameter used was not enough large for the SS phase to 
be observed clearly, as the author of Ref. Xi predicted. 




FIG. 12: (Color online) Ground-state phase diagram at half- 
filling. Circles and squares denote critical points which corre- 
spond to onsets of checkerboard order and superfluid respec- 
tively. The lines are used to guide the eyes. The green region 
between the two lines represents the supersolid (SS) phase. 




0.328 



0.33 0.332 

zt/U 



0.334 



FIG. 13: (Color online) Estimation of quantum critical points 
from intersection of p a L or g for different system sizes. Dashed 
vertical lines are placed at the estimated critical points for the 
CB-SS transition (left) and SS-SF transition (right). 



VI. SUMMARY 



In conclusion, we have investigated the ground-state 
phase diagrams of the 2D extended Bose- Hubbard model 
by performing unbiased quantum Monte Carlo simula- 
tions. Especially, we find that the ground-state phase di- 
agrams by Gutzwiller mean-field approximation are qual- 
itatively correct and the supersolid below and at half- 
filling are stable as well as the 3D system. For the strong 
nearest-neighbor repulsion, we have also confirmed that 
the supersolid phase exits up to large hopping parame- 
ters. Although the 2D result qualitatively agrees with 
the 3D or Gutzwiller mean-field results, the supersolid 
regions shrinks in the lower dimensions due to the quan- 
tum fluctuation. 
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